Регрессия

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(datasets)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:plotly':
## 
##     select

library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 4.1-8

library(ridge)
library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Задание 1

Загрузите данные из файла reglab1.txt. Используя функцию lm, постройте регрессию (используйте разные модели). Выберите наиболее подходящую модель, объясните свой выбор.

data <- read.table("reglab1.txt", header = TRUE)
head(data)

##          z         x         y
## 1 2.836772 0.2710104 0.3083308
## 2 4.987167 0.5895982 0.5149134
## 3 6.412325 0.6517442 0.7304530
## 4 4.641998 0.5819827 0.4614002
## 5 2.793941 0.4636884 0.1911021
## 6 2.934520 0.3835542 0.2612945

n <- dim(data)[1]

data_rand <- data[order(runif(n)),]
df_train <- data_rand[1:as.integer(n*0.8),]
df_test <- data_rand[(as.integer(n*0.8)+1):n,]

plot_ly(data = df_train) %>%
  add_markers(
    z = ~z, x = ~x, y = ~y,
    opacity = 0.7,
    size = I(100)
  ) %>%
  layout(title = "Training Data")

plot_ly(data = df_test) %>%
  add_markers(
    z = ~z, x = ~x, y = ~y,
    opacity = 0.7,
    size = I(100)
  ) %>%
  layout(title = "Test Data")

model1 <- lm(z ~ ., df_train)
predicted_z <- predict(model1, df_test)

model2 <- lm(y ~ ., df_train)
predicted_y <- predict(model2, df_test)

model3 <- lm(x ~ ., df_train)
predicted_x <- predict(model3, df_test)

# Ошибки для моделей
mist_z <- sd(df_test$z - predicted_z)
model1

## 
## Call:
## lm(formula = z ~ ., data = df_train)
## 
## Coefficients:
## (Intercept)            x            y  
##    -0.04412      4.10877      5.01872

cat("Ошибка на тестовых данных: ", mist_z, "\n")

## Ошибка на тестовых данных:  0.3443943

mist_y <- sd(df_test$y - predicted_y)
model2

## 
## Call:
## lm(formula = y ~ ., data = df_train)
## 
## Coefficients:
## (Intercept)            z            x  
##     0.02944      0.18994     -0.77650

cat("Ошибка на тестовых данных: ", mist_y, "\n")

## Ошибка на тестовых данных:  0.07083249

mist_x <- sd(df_test$x - predicted_x)
model3

## 
## Call:
## lm(formula = x ~ ., data = df_train)
## 
## Coefficients:
## (Intercept)            z            y  
##     0.04961      0.22330     -1.11505

cat("Ошибка на тестовых данных: ", mist_x, "\n")

## Ошибка на тестовых данных:  0.08171881

Вывод

Наиболее подходящая модель - модель y, так как меньше ошибка на тестовых данных.

Задание 2

Реализуйте следующий алгоритм для уменьшения количества признаков, используемых для построения регрессии: для каждого k из {0,1,…,d} выбрать подмножество признаков мощности k^1, минимизирующее остаточную сумму квадратов RSS. Используя полученный алгоритм, выберите оптимальное подможество признаков для данных из файла reglab2.txt. Объясните свой выбор. Для генерации всех возможных сочетаний по m элементов из некоторого множества x можно использовать функцию combn(x, m, …).

data <- read.table("reglab2.txt", header = TRUE)
head(data)

##           y         x1         x2        x3         x4
## 1 3.4697200 0.23362753 0.83554899 0.1029653 0.45742768
## 2 0.7684485 0.11791994 0.09054379 0.2587784 0.28395097
## 3 2.8803737 0.09152018 0.79759245 0.1985277 0.69928681
## 4 3.7454851 0.87672239 0.06293480 0.6154150 0.17605768
## 5 1.8539662 0.20740592 0.30348994 0.7759667 0.66735068
## 6 1.2952686 0.04407584 0.36008230 0.6572325 0.04344568

compute_RSS <- function(features, data) {
  formula <- as.formula(paste("y ~", paste(features, collapse = " + ")))
  model <- lm(formula, data = data)
  RSS <- sum(model$residuals^2)
  return(RSS)
}

find_optimal_subset_lm <- function(data) {
  p <- ncol(data) - 1
  features <- colnames(data)[-1]

  best_RSS <- Inf
  best_subset <- NULL

  all_results <- list()

  # Итерация по всем возможным размерам подмножества
  for (m in 1:p) {
    
    subsets <- combn(features, m, simplify = TRUE)

    # Итерация по всем возможным сочетаниям признаков данного размера
    for (i in 1:ncol(subsets)) {
      current_subset <- subsets[, i]
      current_RSS <- compute_RSS(current_subset, data)

      all_results[[length(all_results) + 1]] <- list(subset = current_subset, RSS = current_RSS)

      if (current_RSS < best_RSS) {
        best_RSS <- current_RSS
        best_subset <- current_subset
      }
    }
  }
  
  return(list(all_results = all_results, optimal_subset = best_subset, optimal_RSS = best_RSS))
}

result <- find_optimal_subset_lm(data)

cat("Все подмножества и их RSS:\n")

## Все подмножества и их RSS:

for (i in seq_along(result$all_results)) {
  cat("Подмножество:", paste(result$all_results[[i]]$subset, collapse = ", "), " | RSS:", result$all_results[[i]]$RSS, "\n")
}

## Подмножество: x1  | RSS: 157.2198 
## Подмножество: x2  | RSS: 268.2458 
## Подмножество: x3  | RSS: 393.4905 
## Подмножество: x4  | RSS: 394.5905 
## Подмножество: x1, x2  | RSS: 0.5379617 
## Подмножество: x1, x3  | RSS: 156.3541 
## Подмножество: x1, x4  | RSS: 157.2193 
## Подмножество: x2, x3  | RSS: 267.7955 
## Подмножество: x2, x4  | RSS: 267.8061 
## Подмножество: x3, x4  | RSS: 393.4587 
## Подмножество: x1, x2, x3  | RSS: 0.3322662 
## Подмножество: x1, x2, x4  | RSS: 0.3619682 
## Подмножество: x1, x3, x4  | RSS: 156.3483 
## Подмножество: x2, x3, x4  | RSS: 267.4415 
## Подмножество: x1, x2, x3, x4  | RSS: 0.1928635

cat("\nОптимальное подмножество признаков:", paste(result$optimal_subset, collapse = ", "), "\n")

## 
## Оптимальное подмножество признаков: x1, x2, x3, x4

cat("Максимальная остаточная сумма квадратов (RSS):", result$optimal_RSS, "\n")

## Максимальная остаточная сумма квадратов (RSS): 0.1928635

Задание 3

Загрузите данные из файла cygage.txt. Постройте регрессию, выражающую зависимость возраста исследуемых отложений от глубины залегания, используя веса наблюдений. Оцените качество построенной модели.

data <- read.table("cygage.txt", header = TRUE)
head(data)

##   calAge Depth Weight
## 1      0     0    1.0
## 2   3707    77    0.1
## 3   4150    90    0.1
## 4   5350   107    0.1
## 5   4500   168    0.1
## 6   7260   217    0.1

plot(data, col = 'blue')

model <- lm(calAge ~ Depth, data, weights = data$Weight)

summary <- summary(model)
mid_error = mean(summary$residuals^2)
cat("Среднеквадратичная ошибка (MSE):", mid_error, "\n")

## Среднеквадратичная ошибка (MSE): 227161.2

Задание 4

Загрузите данные Longley (макроэкономические данные). Данные состоят из 7 экономических переменных, наблюдаемых с 1947 по 1962 годы (n=16): GNP.deflator - дефлятор цен, GNP - валовой национальный продукт, Unemployed – число безработных Armed.Forces – число людей в армии Population – население, возраст которого старше 14 лет Year - год Employed – количество занятых Построить регрессию lm(Employed ~ .). Исключите из набора данных longley переменную “Population”. Разделите данные на тестовую и обучающую выборки равных размеров случайным образом. Постройте гребневую регрессию для значений lambda=10^(-3+0.2*i), i=0,…,25, подсчитайте ошибку на тестовой и обучающей выборке для данных значений λ, постройте графики. Объясните полученные результаты.

data(longley)
head(longley)

##      GNP.deflator     GNP Unemployed Armed.Forces Population Year Employed
## 1947         83.0 234.289      235.6        159.0    107.608 1947   60.323
## 1948         88.5 259.426      232.5        145.6    108.632 1948   61.122
## 1949         88.2 258.054      368.2        161.6    109.773 1949   60.171
## 1950         89.5 284.599      335.1        165.0    110.929 1950   61.187
## 1951         96.2 328.975      209.9        309.9    112.075 1951   63.221
## 1952         98.1 346.999      193.2        359.4    113.270 1952   63.639

model <- lm(Employed ~ ., longley)

summary <- summary(model)
mid_error = mean(summary$residuals^2)
cat("Среднеквадратичная ошибка (MSE):", mid_error, "\n")

## Среднеквадратичная ошибка (MSE): 0.0522765

longley$Population <- NULL
n <- dim(longley)[1]

longley_rand <- longley[order(runif(n)),]
df_train <- longley_rand[1:as.integer(n*0.8),]
df_test <- longley_rand[(as.integer(n*0.8)+1):n,]

lambda_values <- 10^seq(-3, 2, by = 0.2)
train_error <- numeric(length(lambda_values))
test_error <- numeric(length(lambda_values))

for (i in seq_along(lambda_values)) {
  lambda <- lambda_values[i]
  model <- lm.ridge(Employed ~ ., df_train, lambda = lambda)

  coef_values <- as.matrix(coef(model))

  train_pred <- as.matrix(cbind(1, df_train[, -ncol(df_train)])) %*% coef_values
  test_pred <- as.matrix(cbind(1, df_test[, -ncol(df_test)])) %*% coef_values

  train_error[i] <- mean((train_pred - df_train$Employed)^2)
  test_error[i] <- mean((test_pred - df_test$Employed)^2)
}

plot(log10(lambda_values), train_error, type = "l", col = "blue", xlab = "log10(lambda)", ylab = "MSE", main = "Гребневая регрессия")
lines(log10(lambda_values), test_error, type = "l", col = "red")
legend("topright", legend = c("Обучающая выборка", "Тестовая выборка"), col = c("blue", "red"), lty = 1)

Из графика видно, что среднеквадратичная ошибка на обучающей выборке больше, чем на тестовой. Причем при увеличении lambda, ошибка растет.

Задание 5

Загрузите данные EuStockMarkets из пакета « datasets». Данные содержат ежедневные котировки на момент закрытия фондовых бирж: Germany DAX (Ibis), Switzerland SMI, France CAC, и UK FTSE. Постройте на одном графике все кривые изменения котировок во времени. Постройте линейную регрессию для каждой модели в отдельности и для всех моделей вместе. Оцените, какая из бирж имеет наибольшую динамику.

data(EuStockMarkets)
head(EuStockMarkets)

##          DAX    SMI    CAC   FTSE
## [1,] 1628.75 1678.1 1772.8 2443.6
## [2,] 1613.63 1688.5 1750.5 2460.2
## [3,] 1606.51 1678.6 1718.0 2448.2
## [4,] 1621.04 1684.1 1708.1 2470.4
## [5,] 1618.16 1686.6 1723.1 2484.7
## [6,] 1610.61 1671.6 1714.3 2466.8

plot(EuStockMarkets[,1], type = "l", xlab="Time", ylab="Closing Price", col="black")
lines(EuStockMarkets[,2], type = "l", col="red")
lines(EuStockMarkets[,3], type = "l", col="green")
lines(EuStockMarkets[,4], type = "l", col="blue")
legend("topleft", legend = colnames(EuStockMarkets), col = 1:4, lty = 1)

model1 <- lm(DAX ~ time(EuStockMarkets), EuStockMarkets)
model2 <- lm(SMI ~ time(EuStockMarkets), EuStockMarkets)
model3 <- lm(CAC ~ time(EuStockMarkets), EuStockMarkets)
model4 <- lm(FTSE ~ time(EuStockMarkets), EuStockMarkets)
model5 <- lm(DAX+SMI+CAC+FTSE ~ time(EuStockMarkets), EuStockMarkets)

cat("-----------Регрессия для DAX----------- \n")

## -----------Регрессия для DAX-----------

model1

## 
## Call:
## lm(formula = DAX ~ time(EuStockMarkets), data = EuStockMarkets)
## 
## Coefficients:
##          (Intercept)  time(EuStockMarkets)  
##            -894557.9                 449.7

cat("-----------Регрессия для SMI----------- \n")

## -----------Регрессия для SMI-----------

model2

## 
## Call:
## lm(formula = SMI ~ time(EuStockMarkets), data = EuStockMarkets)
## 
## Coefficients:
##          (Intercept)  time(EuStockMarkets)  
##           -1428160.2                 717.5

cat("-----------Регрессия для CAC----------- \n")

## -----------Регрессия для CAC-----------

model3

## 
## Call:
## lm(formula = CAC ~ time(EuStockMarkets), data = EuStockMarkets)
## 
## Coefficients:
##          (Intercept)  time(EuStockMarkets)  
##            -405915.3                 204.6

cat("-----------Регрессия для FTSE----------- \n")

## -----------Регрессия для FTSE-----------

model4

## 
## Call:
## lm(formula = FTSE ~ time(EuStockMarkets), data = EuStockMarkets)
## 
## Coefficients:
##          (Intercept)  time(EuStockMarkets)  
##            -865200.4                 435.5

cat("-----------Регрессия для всех----------- \n")

## -----------Регрессия для всех-----------

model5

## 
## Call:
## lm(formula = DAX + SMI + CAC + FTSE ~ time(EuStockMarkets), data = EuStockMarkets)
## 
## Coefficients:
##          (Intercept)  time(EuStockMarkets)  
##             -3593834                  1807

Из результатов видим, что наибольшую динамику имеет SMI, так как у неё самый большой коэффициент при параметре time.

Задание 6

Загрузите данные JohnsonJohnson из пакета «datasets». Данные содержат поквартальную прибыль компании Johnson & Johnson с 1960 по 1980 гг. Постройте на одном графике все кривые изменения прибыли во времени. Постройте линейную регрессию для каждого квартала в отдельности и для всех кварталов вместе. Оцените, в каком квартале компания имеет наибольшую и наименьшую динамику доходности. Сделайте прогноз по прибыли в 2016 году во всех кварталах и в среднем по году.

data(JohnsonJohnson)
JohnsonJohnson

##       Qtr1  Qtr2  Qtr3  Qtr4
## 1960  0.71  0.63  0.85  0.44
## 1961  0.61  0.69  0.92  0.55
## 1962  0.72  0.77  0.92  0.60
## 1963  0.83  0.80  1.00  0.77
## 1964  0.92  1.00  1.24  1.00
## 1965  1.16  1.30  1.45  1.25
## 1966  1.26  1.38  1.86  1.56
## 1967  1.53  1.59  1.83  1.86
## 1968  1.53  2.07  2.34  2.25
## 1969  2.16  2.43  2.70  2.25
## 1970  2.79  3.42  3.69  3.60
## 1971  3.60  4.32  4.32  4.05
## 1972  4.86  5.04  5.04  4.41
## 1973  5.58  5.85  6.57  5.31
## 1974  6.03  6.39  6.93  5.85
## 1975  6.93  7.74  7.83  6.12
## 1976  7.74  8.91  8.28  6.84
## 1977  9.54 10.26  9.54  8.73
## 1978 11.88 12.06 12.15  8.91
## 1979 14.04 12.96 14.85  9.99
## 1980 16.20 14.67 16.02 11.61

time_labels <- time(JohnsonJohnson)[seq(from = 1, to = length(JohnsonJohnson), by = 4)]
col_names <- c('Qtr1','Qtr2','Qtr3','Qtr4')
Qtr1 <- JohnsonJohnson[seq(from = 1, to = length(JohnsonJohnson), by = 4)]
Qtr2 <- JohnsonJohnson[seq(from = 2, to = length(JohnsonJohnson), by = 4)]
Qtr3 <- JohnsonJohnson[seq(from = 3, to = length(JohnsonJohnson), by = 4)]
Qtr4 <- JohnsonJohnson[seq(from = 4, to = length(JohnsonJohnson), by = 4)]

plot(time_labels, Qtr1, type = "l", xlab = "Время", ylab = "Кривая", col = "black")
lines(time_labels, Qtr2, type = "l", col = "red")
lines(time_labels, Qtr3, type = "l", col = "green")
lines(time_labels, Qtr4, type = "l", col = "blue")

legend("topleft", legend = col_names, col = 1:4, lty = 1)

model1 <- lm(Qtr1 ~ time_labels, JohnsonJohnson)
model2 <- lm(Qtr2 ~ time_labels, JohnsonJohnson)
model3 <- lm(Qtr3 ~ time_labels, JohnsonJohnson)
model4 <- lm(Qtr4 ~ time_labels, JohnsonJohnson)
model5 <- lm(Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time_labels, JohnsonJohnson)

cat("-----------Регрессия для Qtr1----------- \n")

## -----------Регрессия для Qtr1-----------

model1

## 
## Call:
## lm(formula = Qtr1 ~ time_labels, data = JohnsonJohnson)
## 
## Coefficients:
## (Intercept)  time_labels  
##   -1364.282        0.695

cat("-----------Регрессия для Qtr2----------- \n")

## -----------Регрессия для Qtr2-----------

model2

## 
## Call:
## lm(formula = Qtr2 ~ time_labels, data = JohnsonJohnson)
## 
## Coefficients:
## (Intercept)  time_labels  
##  -1345.0727       0.6853

cat("-----------Регрессия для Qtr3----------- \n")

## -----------Регрессия для Qtr3-----------

model3

## 
## Call:
## lm(formula = Qtr3 ~ time_labels, data = JohnsonJohnson)
## 
## Coefficients:
## (Intercept)  time_labels  
##  -1382.3170       0.7044

cat("-----------Регрессия для Qtr4----------- \n")

## -----------Регрессия для Qtr4-----------

model4

## 
## Call:
## lm(formula = Qtr4 ~ time_labels, data = JohnsonJohnson)
## 
## Coefficients:
## (Intercept)  time_labels  
##  -1049.5828       0.5349

cat("-----------Регрессия для всех----------- \n")

## -----------Регрессия для всех-----------

model5

## 
## Call:
## lm(formula = Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time_labels, data = JohnsonJohnson)
## 
## Coefficients:
## (Intercept)  time_labels  
##    -5141.25         2.62

Из результатов видим, что наибольшая динамика наблюдается в 3 квартале, так как у него самый большой коэффициент при параметре time. Наименьшая же динамика наблюдается в 4 квартале.

predict1 = coef(model1)[1]+coef(model1)[2]*2016
predict2 = coef(model2)[1]+coef(model2)[2]*2016
predict3 = coef(model3)[1]+coef(model3)[2]*2016
predict4 = coef(model4)[1]+coef(model4)[2]*2016
predict5 = (coef(model5)[1]+coef(model5)[2]*2016)/4

cat("Прогноз прибыли в 2016 году для Qtr1: ", predict1, "\n")

## Прогноз прибыли в 2016 году для Qtr1:  36.75964

cat("Прогноз прибыли в 2016 году для Qtr2: ", predict2, "\n")

## Прогноз прибыли в 2016 году для Qtr2:  36.48945

cat("Прогноз прибыли в 2016 году для Qtr3: ", predict3, "\n")

## Прогноз прибыли в 2016 году для Qtr3:  37.65394

cat("Прогноз прибыли в 2016 году для Qtr4: ", predict4, "\n")

## Прогноз прибыли в 2016 году для Qtr4:  28.79391

cat("Прогноз прибыли в 2016 году в среднем по году: ", predict5, "\n")

## Прогноз прибыли в 2016 году в среднем по году:  34.92424

Задание 7

Загрузите данные sunspot.year из пакета «datasets». Данные содержат количество солнечных пятен с 1700 по 1988 гг. Постройте на графике кривую изменения числа солнечных пятен во времени. Постройте линейную регрессию для данных.

data(sunspot.year)
sunspot.year

## Time Series:
## Start = 1700 
## End = 1988 
## Frequency = 1 
##   [1]   5.0  11.0  16.0  23.0  36.0  58.0  29.0  20.0  10.0   8.0   3.0   0.0
##  [13]   0.0   2.0  11.0  27.0  47.0  63.0  60.0  39.0  28.0  26.0  22.0  11.0
##  [25]  21.0  40.0  78.0 122.0 103.0  73.0  47.0  35.0  11.0   5.0  16.0  34.0
##  [37]  70.0  81.0 111.0 101.0  73.0  40.0  20.0  16.0   5.0  11.0  22.0  40.0
##  [49]  60.0  80.9  83.4  47.7  47.8  30.7  12.2   9.6  10.2  32.4  47.6  54.0
##  [61]  62.9  85.9  61.2  45.1  36.4  20.9  11.4  37.8  69.8 106.1 100.8  81.6
##  [73]  66.5  34.8  30.6   7.0  19.8  92.5 154.4 125.9  84.8  68.1  38.5  22.8
##  [85]  10.2  24.1  82.9 132.0 130.9 118.1  89.9  66.6  60.0  46.9  41.0  21.3
##  [97]  16.0   6.4   4.1   6.8  14.5  34.0  45.0  43.1  47.5  42.2  28.1  10.1
## [109]   8.1   2.5   0.0   1.4   5.0  12.2  13.9  35.4  45.8  41.1  30.1  23.9
## [121]  15.6   6.6   4.0   1.8   8.5  16.6  36.3  49.6  64.2  67.0  70.9  47.8
## [133]  27.5   8.5  13.2  56.9 121.5 138.3 103.2  85.7  64.6  36.7  24.2  10.7
## [145]  15.0  40.1  61.5  98.5 124.7  96.3  66.6  64.5  54.1  39.0  20.6   6.7
## [157]   4.3  22.7  54.8  93.8  95.8  77.2  59.1  44.0  47.0  30.5  16.3   7.3
## [169]  37.6  74.0 139.0 111.2 101.6  66.2  44.7  17.0  11.3  12.4   3.4   6.0
## [181]  32.3  54.3  59.7  63.7  63.5  52.2  25.4  13.1   6.8   6.3   7.1  35.6
## [193]  73.0  85.1  78.0  64.0  41.8  26.2  26.7  12.1   9.5   2.7   5.0  24.4
## [205]  42.0  63.5  53.8  62.0  48.5  43.9  18.6   5.7   3.6   1.4   9.6  47.4
## [217]  57.1 103.9  80.6  63.6  37.6  26.1  14.2   5.8  16.7  44.3  63.9  69.0
## [229]  77.8  64.9  35.7  21.2  11.1   5.7   8.7  36.1  79.7 114.4 109.6  88.8
## [241]  67.8  47.5  30.6  16.3   9.6  33.2  92.6 151.6 136.3 134.7  83.9  69.4
## [253]  31.5  13.9   4.4  38.0 141.7 190.2 184.8 159.0 112.3  53.9  37.5  27.9
## [265]  10.2  15.1  47.0  93.8 105.9 105.5 104.5  66.6  68.9  38.0  34.5  15.5
## [277]  12.6  27.5  92.5 155.4 154.7 140.5 115.9  66.6  45.9  17.9  13.4  29.2
## [289] 100.2

plot(sunspot.year, col='blue')

model <- lm(sunspot.year[seq(1,length(sunspot.year),1)] ~ seq(1700,1988,1), sunspot.year)

cat("-----------Регрессия----------- \n")

## -----------Регрессия-----------

model

## 
## Call:
## lm(formula = sunspot.year[seq(1, length(sunspot.year), 1)] ~ 
##     seq(1700, 1988, 1), data = sunspot.year)
## 
## Coefficients:
##        (Intercept)  seq(1700, 1988, 1)  
##         -130.42119             0.09709

summary <- summary(model)
mid_error = mean(summary$residuals^2)
cat("Среднеквадратичная ошибка (MSE):", mid_error, "\n")

## Среднеквадратичная ошибка (MSE): 1487.204

Задание 8

Загрузите данные из файла пакета «UKgas.scv». Данные содержат объемы ежеквартально потребляемого газа в Великобритании с 1960 по 1986 гг. Постройте линейную регрессию для каждого квартала в отдельности и для всех кварталов вместе. Оцените, в каком квартале потребление газа имеет наибольшую и наименьшую динамику доходности. Сделайте прогноз по потреблению газа в 2016 году во всех кварталах и в среднем по году.

data <- read.csv("UKgas.csv", stringsAsFactors = TRUE)
head(data)

##   X    time UKgas
## 1 1 1960.00 160.1
## 2 2 1960.25 129.7
## 3 3 1960.50  84.8
## 4 4 1960.75 120.1
## 5 5 1961.00 160.1
## 6 6 1961.25 124.9

Qtr1 <- UKgas[seq(from = 1, to = dim(data)[1], by = 4)]
Qtr2 <- UKgas[seq(from = 2, to = dim(data)[1], by = 4)]
Qtr3 <- UKgas[seq(from = 3, to = dim(data)[1], by = 4)]
Qtr4 <- UKgas[seq(from = 4, to = dim(data)[1], by = 4)]

model1 <- lm(Qtr1 ~ time[seq(from = 1, to = dim(data)[1], by = 4)], data)
model2 <- lm(Qtr2 ~ time[seq(from = 2, to = dim(data)[1], by = 4)], data)
model3 <- lm(Qtr3 ~ time[seq(from = 3, to = dim(data)[1], by = 4)], data)
model4 <- lm(Qtr4 ~ time[seq(from = 4, to = dim(data)[1], by = 4)], data)
model5 <- lm(Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time[seq(from = 1, to = dim(data)[1], by = 4)], data)

cat("-----------Регрессия для Qtr1----------- \n")

## -----------Регрессия для Qtr1-----------

model1

## 
## Call:
## lm(formula = Qtr1 ~ time[seq(from = 1, to = dim(data)[1], by = 4)], 
##     data = data)
## 
## Coefficients:
##                                    (Intercept)  
##                                      -78854.23  
## time[seq(from = 1, to = dim(data)[1], by = 4)]  
##                                          40.22

cat("-----------Регрессия для Qtr2----------- \n")

## -----------Регрессия для Qtr2-----------

model2

## 
## Call:
## lm(formula = Qtr2 ~ time[seq(from = 2, to = dim(data)[1], by = 4)], 
##     data = data)
## 
## Coefficients:
##                                    (Intercept)  
##                                      -35297.23  
## time[seq(from = 2, to = dim(data)[1], by = 4)]  
##                                          18.04

cat("-----------Регрессия для Qtr3----------- \n")

## -----------Регрессия для Qtr3-----------

model3

## 
## Call:
## lm(formula = Qtr3 ~ time[seq(from = 3, to = dim(data)[1], by = 4)], 
##     data = data)
## 
## Coefficients:
##                                    (Intercept)  
##                                      -15403.73  
## time[seq(from = 3, to = dim(data)[1], by = 4)]  
##                                           7.89

cat("-----------Регрессия для Qtr4----------- \n")

## -----------Регрессия для Qtr4-----------

model4

## 
## Call:
## lm(formula = Qtr4 ~ time[seq(from = 4, to = dim(data)[1], by = 4)], 
##     data = data)
## 
## Coefficients:
##                                    (Intercept)  
##                                      -59112.72  
## time[seq(from = 4, to = dim(data)[1], by = 4)]  
##                                          30.14

cat("-----------Регрессия для всех----------- \n")

## -----------Регрессия для всех-----------

model5

## 
## Call:
## lm(formula = Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time[seq(from = 1, to = dim(data)[1], 
##     by = 4)], data = data)
## 
## Coefficients:
##                                    (Intercept)  
##                                     -188636.85  
## time[seq(from = 1, to = dim(data)[1], by = 4)]  
##                                          96.29

Из результатов видим, что наибольшая динамика наблюдается в 1 квартале, так как у него самый большой коэффициент при параметре time. Наименьшая же динамика наблюдается в 3 квартале.

predict1 = coef(model1)[1]+coef(model1)[2]*2016
predict2 = coef(model2)[1]+coef(model2)[2]*2016
predict3 = coef(model3)[1]+coef(model3)[2]*2016
predict4 = coef(model4)[1]+coef(model4)[2]*2016
predict5 = (coef(model5)[1]+coef(model5)[2]*2016)/4

cat("Прогноз прибыли в 2016 году для Qtr1: ", predict1, "\n")

## Прогноз прибыли в 2016 году для Qtr1:  2230.936

cat("Прогноз прибыли в 2016 году для Qtr2: ", predict2, "\n")

## Прогноз прибыли в 2016 году для Qtr2:  1072.375

cat("Прогноз прибыли в 2016 году для Qtr3: ", predict3, "\n")

## Прогноз прибыли в 2016 году для Qtr3:  501.9919

cat("Прогноз прибыли в 2016 году для Qtr4: ", predict4, "\n")

## Прогноз прибыли в 2016 году для Qtr4:  1654.785

cat("Прогноз прибыли в 2016 году в среднем по году: ", predict5, "\n")

## Прогноз прибыли в 2016 году в среднем по году:  1372.787

Задание 9

Загрузите данные cars из пакета «datasets». Данные содержат зависимости тормозного пути автомобиля (футы) от его скорости (мили в час). Данные получены в 1920 г. Постройте регрессионную модель и оцените длину тормозного пути при скорости 40 миль в час.

data(cars)
head(cars)

##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10

model <- lm(dist ~ speed, cars)

cat("-----------Регрессия----------- \n")

## -----------Регрессия-----------

model

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Coefficients:
## (Intercept)        speed  
##     -17.579        3.932

predict = coef(model)[1]+coef(model)[2]*40
cat("Прогноз длины тормозного пути при скорости 40 миль в час: ", predict, "\n")

## Прогноз длины тормозного пути при скорости 40 миль в час:  139.7173